Airline tickets: customers loyalty

Report tecnico progetto Esame Analisi dei Big Data

Il seguente progetto riguarda l’analisi un dataset di biglietti aerei venduti online da una compagnia aerea che effettua collegamenti dalle principali città italiane alle principali città europee.

Analisi

Il presente progetto contiene i seguenti punti:

Analisi descrittiva del dataset:

• 1. Descrizione dei dati osservati (indici statistici univariati e/o bivariati, rappresentazioni grafiche, ecc.);
• 2. Segmentazione della clientela sulla base delle variabili osservate;
• 3. Suddivisione casuale del dataset in campione di apprendimento (2/3 · n) e campione test (1/3 · n). Stima, sul campione di apprendimento, di 3 modelli (Lasso, PCR, PLS ) in funzione della variabile “Cfidelity” (-> factor, levels:alta, bassa. Creata tramite ricodifica della Variabile “Fidelity”).

Importazione del dataset

load("~/Documents/r Wdirectory/Dati_Caria.RData")
summary(dati)

##       Data             Unit_Ticket     Seat       PriorBoard  Luggage    
##  Min.   :2010-01-01   Min.   : 87.9   No :11765   No :11576   No :13111  
##  1st Qu.:2012-06-21   1st Qu.:158.0   Yes:12233   Yes:12422   Yes:10887  
##  Median :2014-12-27   Median :173.8                                      
##  Mean   :2014-12-28   Mean   :175.3                                      
##  3rd Qu.:2017-06-30   3rd Qu.:191.5                                      
##  Max.   :2019-12-31   Max.   :274.6                                      
##  Return         Fidelity         BookTime         Taxes          Discount     
##  No :12897   Min.   :  0.00   Min.   : 3.00   Min.   : 8.88   Min.   : 0.000  
##  Yes:11101   1st Qu.: 33.35   1st Qu.: 9.00   1st Qu.:20.23   1st Qu.: 0.860  
##              Median : 43.06   Median :11.00   Median :23.33   Median : 1.940  
##              Mean   : 45.31   Mean   :11.76   Mean   :23.68   Mean   : 3.493  
##              3rd Qu.: 56.29   3rd Qu.:15.00   3rd Qu.:27.00   3rd Qu.: 6.210  
##              Max.   :100.00   Max.   :22.00   Max.   :41.32   Max.   :15.160  
##       NPax          Departure         Arrival            ModPag     
##  Min.   : 1.000   Milan  :5481   London   :11212   C_debito : 3095  
##  1st Qu.: 4.000   Rome   :9858   Paris    :  192   C_credito:20903  
##  Median : 6.000   Naples : 380   Madrid   : 2938                    
##  Mean   : 5.726   Palermo:2806   Amsterdam: 9656                    
##  3rd Qu.: 7.000   Bergamo:1656                                      
##  Max.   :12.000   Venice :3817                                      
##    NBookMonth     NAccWebWeek    NComplaintsYear  NRefundYear    
##  Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   : 0.000  
##  1st Qu.:3.000   1st Qu.:2.000   1st Qu.:1.000   1st Qu.: 0.000  
##  Median :4.000   Median :3.000   Median :1.000   Median : 0.000  
##  Mean   :3.993   Mean   :2.996   Mean   :1.157   Mean   : 5.853  
##  3rd Qu.:5.000   3rd Qu.:4.000   3rd Qu.:2.000   3rd Qu.: 0.000  
##  Max.   :9.000   Max.   :7.000   Max.   :6.000   Max.   :65.000  
##   NCancelYear   
##  Min.   : 1.00  
##  1st Qu.: 4.00  
##  Median : 5.00  
##  Mean   : 5.01  
##  3rd Qu.: 6.00  
##  Max.   :10.00

str(dati)

## 'data.frame':    23998 obs. of  19 variables:
##  $ Data           : Date, format: "2010-01-01" "2010-01-01" ...
##  $ Unit_Ticket    : num  174 160 173 189 152 ...
##  $ Seat           : Factor w/ 2 levels "No","Yes": 2 1 2 2 1 1 2 1 1 2 ...
##  $ PriorBoard     : Factor w/ 2 levels "No","Yes": 2 2 1 1 2 1 2 2 2 1 ...
##  $ Luggage        : Factor w/ 2 levels "No","Yes": 2 1 1 1 2 1 1 2 2 2 ...
##  $ Return         : Factor w/ 2 levels "No","Yes": 1 1 2 1 2 1 1 1 1 1 ...
##  $ Fidelity       : num  53.5 47.9 39 40.9 30.1 ...
##  $ BookTime       : num  18 12 10 14 13 15 14 18 18 16 ...
##  $ Taxes          : num  24.6 24.5 18.7 19.6 14.9 ...
##  $ Discount       : num  6.97 1.23 0.7 1.01 0.48 1.33 1.06 9.91 8.97 8.02 ...
##  $ NPax           : num  4 6 4 9 4 5 5 6 7 7 ...
##  $ Departure      : Factor w/ 6 levels "Milan","Rome",..: 2 2 1 6 2 2 2 2 5 4 ...
##  $ Arrival        : Factor w/ 4 levels "London","Paris",..: 1 1 4 1 1 1 1 1 1 1 ...
##  $ ModPag         : Factor w/ 2 levels "C_debito","C_credito": 2 2 2 2 2 2 2 1 2 2 ...
##  $ NBookMonth     : num  5 4 3 5 5 4 5 4 4 4 ...
##  $ NAccWebWeek    : num  3 3 3 5 3 2 4 3 2 3 ...
##  $ NComplaintsYear: num  1 0 2 3 1 1 2 0 0 1 ...
##  $ NRefundYear    : num  0 0 27 0 0 0 28 0 14 0 ...
##  $ NCancelYear    : num  3 4 6 6 4 4 5 5 4 5 ...

sum(is.na(dati)) #non sono presenti valori mancanti

## [1] 0

Per prima cosa, è stato importato il dataset. Dal summary si può osservare che non sono presenti valori mancanti. Dalla funzione str() si osserva che ci sono 11 variabili num e 7 factor.

Analisi esplorativa del dataset: grafici

library("tidyverse")

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.6     ✓ purrr   0.3.4
## ✓ tibble  3.1.4     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

#1
tab <-dati %>% 
  count(Arrival) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(Arrival= reorder(Arrival, proportion))

tab %>%
  ggplot(aes(Arrival, proportion))+ 
  geom_bar(stat = "identity") 

La proporzione di biglietti venduti per le 4 diverse destinazioni

tab2 <-dati %>% 
  count(Departure) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(Departure= reorder(Departure, proportion))


tab2 %>%
  ggplot(aes(Departure, proportion))+ 
  geom_bar(stat = "identity") 

percentuale di biglietti venduti in base ai 6 diversi aereoporti di provenienza

library("gridExtra")

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

tab3 <- dati %>% 
  count(ModPag) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(ModPag = reorder(ModPag, proportion))

tab3<-tab3 %>%         
  ggplot(aes(ModPag, proportion))+
  geom_bar(stat = "identity") # quasi il 90% (0.87) dei clienti ha utilizzato 
#la carta di credito per il pagamento

#4
tab4 <- dati %>% 
  count(Luggage) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(Luggage = reorder(Luggage, proportion))

tab4<- tab4 %>%         
  ggplot(aes(Luggage, proportion))+
  geom_bar(stat = "identity") # il 54% dei clienti ha effettuato una 
#prenotazione senza imbarcare un bagaglio aggiuntivo  
grid.arrange(tab3, tab4, ncol=2)

Quasi il 90% (0.87) dei clienti ha utilizzato la carta di credito per il pagamento, mentre dal grafico a sinistra osserviamo che il 54% dei clienti ha effettuato una prenotazione senza imbarcare un bagaglio aggiuntivo

#5
tab5 <- dati %>% 
  count(PriorBoard) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(PriorBoard = reorder(PriorBoard, proportion))

tab5<- tab5 %>%         
  ggplot(aes(PriorBoard, proportion))+
  geom_bar(stat = "identity") # il 52% (0.517) dei clienti ha effettuato una 
#prenotazione con imbarco prioritario

#6
tab6 <- dati %>% 
  count(Seat) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(Seat = reorder(Seat, proportion))

tab6<- tab6 %>%         
  ggplot(aes(Seat, proportion))+
  geom_bar(stat = "identity") # il 51% (0.509) dei clienti ha effettuato una 
#prenotazione del posto a sedere



#7
tab7 <- dati %>% 
  count(Return) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(Return = reorder(Return, proportion))

tab7<-tab7 %>%         
  ggplot(aes(Return, proportion))+
  geom_bar(stat = "identity") # il 54% (0.537) dei clienti ha effettuato una 
#prenotazione solo andata
grid.arrange(tab5, tab6, tab7, ncol=3)

Il primo grafico a sinistra mostra che il 52% (0.517) dei clienti ha effettuato una prenotazione con imbarco prioritario. Dal secondo grafico si osserva che il 51% (0.509) dei clienti ha effettuato una prenotazione del posto a sedere. Dal terzo grafico invece, si osserva che il 54% (0.537) dei clienti ha effettuato una prenotazione solo andata.

filter(dati) %>% 
  ggplot(aes(Unit_Ticket, Fidelity, color= Arrival)) +
  geom_point() 

Il prezzo dei biglietti di Londra è più alto rispetto a quello di Amsterdam e Madrid. Il costo inferiore sembra essere dovuto al fatto che l’indice di fidelizzazione del cliente è basso per quest’ultime, mentre assume valori alti per Londra.

filter(dati) %>% 
  ggplot(aes(Fidelity, Discount, color= Arrival)) +
  geom_point() 

I clienti che acquistano i biglietti per Londra risultano essere più fedeli rispetto a quelli di Amsterdam e Madrid e per questo godono di uno sconto maggiore a dispetto di quest’ultimi

dati  %>% 
  ggplot(aes( y=Fidelity, Departure)) + 
  geom_boxplot(coef=3) +
  geom_jitter(width= 0.1, alpha =0.2)

l’indice di fedeltà risulta essere mediamente più alto per le prenotazioni con partenza da Venezia

dati  %>% 
  ggplot(aes( y=Fidelity, Arrival)) + 
  geom_boxplot(coef=3) +
  geom_jitter(width= 0.1, alpha =0.2)

Per quanto riguarda l’areoporto di Arrivo, l’indice di fedeltà risulta essere mediamente più alto per le prenotazioni con arrivo a Londra

library("lubridate")

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

dates<-ymd(dati$Data)


month_day_year<-tibble(date = dates, 
                       month = month(dates, label = TRUE),
                       day = day(dates),
                       year = year(dates)) 


head(month_day_year)

#11  graf
dati2= dati 
dati2<- mutate(dati2, month= month_day_year$month ) #inserisco una nuova colonna 
dati2  %>% 
  ggplot(aes(month, Unit_Ticket)) + 
  geom_boxplot(coef=3) +
  geom_jitter(width= 0.1, alpha =0.2)

#le distribuzioni del prezzo del biglietto nel corso dell'anno è molto simile. Il prezzo del biglietto più basso si è registrato a Novembre, il più prezzo più alto è stato registato a maggio, quello più basso a novembre.

Dopo aver creato un oggetto tibble per dividere giorno, mese e anno, inserisco una nuova colonna in una copia del dataset (dati2). tramite ggplot, ottengodei boxplot riguardanti le distribuzioni del prezzo del biglietto nel corso dell’anno. tali distribuzioni risultano essere molto simile. Il prezzo del biglietto più basso si è registrato a Novembre, il più prezzo più alto è stato registato a maggio.

#ricodifica dell'indice di fedeltà
dati2$clientFidelity[dati2$Fidelity<=50] = "basso"
dati2$clientFidelity[dati2$Fidelity>50] = "alto"

fedeltà <-as.factor(dati2$clientFidelity)
class(fedeltà)

## [1] "factor"

dati2<- dati2 %>% mutate(Cfidelity= fedeltà)

#eliminazione colonne rindondanti dopo la ricodifica
dati2= dati2[-21] #elimino ClientFidelity
dati2= dati2[-7] #elimino la vecchia colonna Fidelity
dati2= dati2[-19]



#12
tab12 <- dati2 %>% 
  count(Cfidelity) %>%
  mutate(proportion= n/sum(n)) %>%
  mutate(Cfidelity = reorder(Cfidelity, proportion))

tab12 %>%         
  ggplot(aes(Cfidelity, proportion))+
  geom_bar(stat = "identity") 

é stato ricodificato l’indice di fedelizzazione: è stata creata una nuova variabile Cfidelity) in cui se l’indice di fidelizzazione (tra 0 e 100) è maggiore di 50 verrà assegnato il valore “alto”, se invece sarà minore di 50 “basso”. Il 66% (0.659) dei clienti ha un indice di fidelizzazione basso (minore al 50) mentre il 34 % dei clienti ha un valore alto per questo indice (maggiore di 50)

library(dslabs)


arrival_departures <- count(dati, Departure, Arrival) #creo un nuovo oggetto 
#tibble in cui conto i viaggi in base a aereoporto di partenza e aereoporto 
#di arrivo

arrival_departures

summary(arrival_departures)

##    Departure      Arrival        n         
##  Milan  :4   London   :6   Min.   :   2.0  
##  Rome   :4   Paris    :6   1st Qu.:  81.0  
##  Naples :4   Madrid   :6   Median : 387.5  
##  Palermo:4   Amsterdam:6   Mean   : 999.9  
##  Bergamo:4                 3rd Qu.:1150.8  
##  Venice :4                 Max.   :4448.0

#13
p1<-arrival_departures %>% unite(flights, Departure, Arrival) #unisco le colonne 
#Unite tramite la funzione Unite della libreria Tidyr
p1

p1 %>%
  mutate( flights= reorder(flights, n)) %>%
  ggplot(aes(flights, n)) +
  geom_bar(stat= "identity")+
  coord_flip()+
  theme(axis.text.y = element_text(size = 8)) +
  xlab("")

Si osserva che i biglietti maggiormente venduti nel sito sono per il volo Roma-Londra, seguito da Milano-Amsterdam e Roma-Amsterdam

index<- dati2$Unit_Ticket; x<- dati2$Unit_Ticket[index]
m <-mean(x) ; s<- sd(x)
c(average=m, sd= s )

##   average        sd 
## 177.70689  25.81076

norm_dist<- data.frame(x= seq(-4, 4, len=50)*s + m) %>% 
  mutate(density= dnorm(x, m, s))

#14
dati2 %>% ggplot(aes(Unit_Ticket))+
  geom_histogram(aes(y=..density..), binwidth = 1, color = "black")+
  geom_line(aes(x, density), data= norm_dist, lwd=1.5)

La distribuzione del prezzo dei biglietti venduti sembra essere approssimabile a una normale.

PCA

library(FactoMineR) 

attach(dati)
new_dati<- tibble(Unit_Ticket, Fidelity, BookTime, Taxes, Discount, NPax, NBookMonth,NAccWebWeek, NComplaintsYear, NRefundYear,NCancelYear, Seat, PriorBoard, Luggage, Return, Departure, Arrival, ModPag)

new_dati2 <- tibble(Unit_Ticket, Fidelity, BookTime, Taxes, Discount, NPax, NBookMonth,NAccWebWeek, NComplaintsYear, NRefundYear, NCancelYear, Arrival )
res <- PCA(new_dati2,
           scale.unit = TRUE,
           ind.sup = 12000: 23998,
           quanti.sup = 6:10,
           quali.sup = 12,
           ncp=8)

summary(res)

## 
## Call:
## PCA(X = new_dati2, scale.unit = TRUE, ncp = 8, ind.sup = 12000:23998,  
##      quanti.sup = 6:10, quali.sup = 12) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6
## Variance               3.546   1.000   0.585   0.470   0.229   0.169
## % of var.             59.097  16.671   9.758   7.839   3.811   2.824
## Cumulative % of var.  59.097  75.768  85.525  93.365  97.176 100.000
## 
## Individuals (the 10 first)
##                     Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1               |  2.690 |  1.566  0.006  0.339 | -1.739  0.025  0.418 | -0.612
## 2               |  1.273 | -0.372  0.000  0.085 | -0.864  0.006  0.461 | -0.264
## 3               |  1.685 | -1.231  0.004  0.534 |  0.840  0.006  0.248 |  0.731
## 4               |  1.658 | -0.261  0.000  0.025 |  0.837  0.006  0.255 |  0.952
## 5               |  2.605 | -1.827  0.008  0.492 | -0.898  0.007  0.119 |  0.214
## 6               |  2.499 |  1.208  0.003  0.234 | -0.902  0.007  0.130 |  1.827
## 7               |  1.666 | -0.867  0.002  0.271 | -0.024  0.000  0.000 |  0.173
## 8               |  3.457 |  3.270  0.025  0.895 | -0.016  0.000  0.000 |  0.395
## 9               |  3.018 |  2.754  0.018  0.833 | -0.881  0.006  0.085 |  0.127
## 10              |  1.935 |  1.632  0.006  0.711 | -0.017  0.000  0.000 | -0.133
##                    ctr   cos2  
## 1                0.005  0.052 |
## 2                0.001  0.043 |
## 3                0.008  0.188 |
## 4                0.013  0.330 |
## 5                0.001  0.007 |
## 6                0.047  0.534 |
## 7                0.000  0.011 |
## 8                0.002  0.013 |
## 9                0.000  0.002 |
## 10               0.000  0.005 |
## 
## Supplementary individuals (the 10 first)
##                     Dist    Dim.1   cos2    Dim.2   cos2    Dim.3   cos2  
## 12000           |  1.403 | -0.515  0.135 |  0.006  0.000 | -0.310  0.049 |
## 12001           |  1.325 | -0.638  0.232 | -0.028  0.000 |  0.647  0.239 |
## 12002           |  0.913 | -0.207  0.051 | -0.014  0.000 |  0.331  0.131 |
## 12003           |  2.435 | -2.004  0.678 | -0.017  0.000 | -0.388  0.025 |
## 12004           |  2.964 |  2.523  0.725 |  0.891  0.090 | -1.105  0.139 |
## 12005           |  1.726 | -1.112  0.415 | -0.853  0.244 | -0.672  0.152 |
## 12006           |  1.462 | -0.825  0.318 |  0.847  0.336 |  0.441  0.091 |
## 12007           |  2.831 |  2.742  0.938 | -0.003  0.000 | -0.158  0.003 |
## 12008           |  2.336 | -1.759  0.567 | -0.898  0.148 |  0.698  0.089 |
## 12009           |  3.771 |  3.465  0.844 |  0.855  0.051 |  1.005  0.071 |
## 
## Variables
##                    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Unit_Ticket     |  0.793 17.726  0.629 | -0.010  0.009  0.000 |  0.544 50.632
## Fidelity        |  0.902 22.948  0.814 | -0.005  0.002  0.000 |  0.223  8.479
## BookTime        |  0.889 22.283  0.790 |  0.001  0.000  0.000 | -0.109  2.041
## Taxes           |  0.765 16.490  0.585 |  0.022  0.047  0.000 | -0.382 24.968
## Discount        |  0.854 20.553  0.729 | -0.001  0.000  0.000 | -0.285 13.845
## NCancelYear     | -0.004  0.000  0.000 |  1.000 99.942  1.000 |  0.014  0.034
##                   cos2  
## Unit_Ticket      0.296 |
## Fidelity         0.050 |
## BookTime         0.012 |
## Taxes            0.146 |
## Discount         0.081 |
## NCancelYear      0.000 |
## 
## Supplementary continuous variables
##                    Dim.1   cos2    Dim.2   cos2    Dim.3   cos2  
## NPax            |  0.558  0.311 | -0.002  0.000 |  0.198  0.039 |
## NBookMonth      |  0.001  0.000 | -0.483  0.233 | -0.018  0.000 |
## NAccWebWeek     | -0.006  0.000 |  0.358  0.128 |  0.004  0.000 |
## NComplaintsYear | -0.009  0.000 |  0.357  0.128 |  0.001  0.000 |
## NRefundYear     |  0.003  0.000 |  0.017  0.000 |  0.007  0.000 |
## 
## Supplementary categories
##                      Dist     Dim.1    cos2  v.test     Dim.2    cos2  v.test  
## London          |   1.429 |   1.409   0.973  76.366 |  -0.008   0.000  -0.846 |
## Paris           |   0.730 |  -0.567   0.605  -2.949 |  -0.038   0.003  -0.376 |
## Madrid          |   0.796 |  -0.715   0.807 -15.489 |   0.025   0.001   1.022 |
## Amsterdam       |   1.421 |  -1.388   0.954 -66.715 |   0.003   0.000   0.247 |
##                   Dim.3    cos2  v.test  
## London           -0.027   0.000  -3.548 |
## Paris             0.302   0.171   3.864 |
## Madrid            0.246   0.096  13.140 |
## Amsterdam        -0.049   0.001  -5.850 |

La PCA può essere utile per la segmentazione della clientela. Si osserva che il primo asse fattoriale è correlato positivamente con le variabili : Unit_Tickets, Fidelity, BookTime, Taxes, Discount. Il secondo asse fattoriale è correlato positivamente con NCancelYear. inoltre sono state proiettate le variabili illustrative.

plot(res, cex=0.8, habillage="Arrival", label = "none")

Si osserva che i punti che i punti sul primo quadrante riguardano la tipologia di clienti con un alto indice di fidelizzazione per le quali prenotazioni si caratterizzano per il costo del biglietto elevato, per uno sconto alto sulla prenotazione, e per un alto numero di giorni di anticipo sulla prenotazione rispetto alla partenza.

Suddivisione casuale del dataset in training set e test set

set.seed(123)
x=model.matrix(Cfidelity ~.,dati2)[,-1]

train = sample(nrow(x), nrow(x)*2/3)
test=(-train)
dati_test = dati2[-train, ] #creo il test set
dati_train = dati2[train,]

y=as.numeric(dati2$Cfidelity)
y.test=as.numeric(y[test])

Il dataset è stato diviso casualmente per 2/3 in train e per 1/3 test set

#lasso

library(glmnet)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loaded glmnet 4.0-2

grid=10^seq(10,-2,length=100)
lasso.mod=glmnet(x[train ,],as.numeric(y[train]),alpha=1,lambda=grid) 
plot(lasso.mod)

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

cv.out= cv.glmnet(x[train ,],as.numeric(y[train]),alpha=1)
plot(cv.out)

bestlam=cv.out$lambda.min
lasso.pred=predict(lasso.mod,s=bestlam ,newx=x[test,])
table(y.test, round(lasso.pred,0)) 

##       
## y.test    1    2
##      1 1984  702
##      2   76 5238

778 malclassificati su 8000 –> 0.0971 –>9,7%

out=glmnet(x,y,alpha=1,lambda=grid)
lasso.coef=predict(out,type="coefficients",s=bestlam)[1:19,] 
lasso.coef 

##      (Intercept)             Data      Unit_Ticket          SeatYes 
##      2.573904903      0.000000000     -0.003758809      0.001333084 
##    PriorBoardYes       LuggageYes        ReturnYes         BookTime 
##      0.000000000     -0.100771968      0.246232407     -0.006281931 
##            Taxes         Discount             NPax    DepartureRome 
##     -0.004634191     -0.044389692      0.000000000      0.000000000 
##  DepartureNaples DeparturePalermo DepartureBergamo  DepartureVenice 
##      0.000000000      0.000000000      0.000000000      0.000000000 
##     ArrivalParis    ArrivalMadrid ArrivalAmsterdam 
##      0.000000000      0.010621467      0.000000000

lasso.coef[lasso.coef!=0] 

##   (Intercept)   Unit_Ticket       SeatYes    LuggageYes     ReturnYes 
##   2.573904903  -0.003758809   0.001333084  -0.100771968   0.246232407 
##      BookTime         Taxes      Discount ArrivalMadrid 
##  -0.006281931  -0.004634191  -0.044389692   0.010621467

solamente 8 dei coefficenti stimati hanno valore diverso da 0. Per cui il lasso model scelto dalla cross-validation coniene solo 8 variabili

summary(lasso.mod)

##           Length Class     Mode   
## a0         100   -none-    numeric
## beta      2400   dgCMatrix S4     
## df         100   -none-    numeric
## dim          2   -none-    numeric
## lambda     100   -none-    numeric
## dev.ratio  100   -none-    numeric
## nulldev      1   -none-    numeric
## npasses      1   -none-    numeric
## jerr         1   -none-    numeric
## offset       1   -none-    logical
## call         5   -none-    call   
## nobs         1   -none-    numeric

#PCR

library(pls)

## 
## Attaching package: 'pls'

## The following object is masked from 'package:stats':
## 
##     loadings

pcr_model <- pcr(as.numeric(Cfidelity) ~ ., data = dati2, 
                 scale = TRUE, 
                 validation = "CV")

summary(pcr_model)

## Data:    X dimension: 23998 24 
##  Y dimension: 23998 1
## Fit method: svdpc
## Number of components considered: 24
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV          0.4738   0.3038   0.3038   0.3038   0.3038   0.2961   0.2956
## adjCV       0.4738   0.3038   0.3038   0.3038   0.3038   0.2961   0.2956
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## CV      0.2925   0.2924   0.2921    0.2921    0.2919    0.2917    0.2916
## adjCV   0.2925   0.2924   0.2921    0.2921    0.2919    0.2916    0.2916
##        14 comps  15 comps  16 comps  17 comps  18 comps  19 comps  20 comps
## CV       0.2916    0.2916    0.2916    0.2913    0.2867    0.2866    0.2809
## adjCV    0.2916    0.2916    0.2916    0.2913    0.2867    0.2866    0.2809
##        21 comps  22 comps  23 comps  24 comps
## CV         0.28      0.28    0.2792    0.2792
## adjCV      0.28      0.28    0.2792    0.2792
## 
## TRAINING: % variance explained
##                        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## X                        22.04    28.77    35.01    40.94    46.53    51.43
## as.numeric(Cfidelity)    58.90    58.90    58.91    58.92    60.96    61.08
##                        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
## X                        56.24    60.79    65.09     69.34     73.55     77.69
## as.numeric(Cfidelity)    61.91    61.95    62.03     62.03     62.10     62.16
##                        13 comps  14 comps  15 comps  16 comps  17 comps
## X                         81.82     85.38     88.70     91.69     94.06
## as.numeric(Cfidelity)     62.17     62.17     62.17     62.17     62.26
##                        18 comps  19 comps  20 comps  21 comps  22 comps
## X                         95.85     97.45     98.55     99.15     99.56
## as.numeric(Cfidelity)     63.45     63.46     64.90     65.13     65.13
##                        23 comps  24 comps
## X                         99.81    100.00
## as.numeric(Cfidelity)     65.34     65.34

# Plot the cross validation MSE
validationplot(pcr_model, val.type="MSEP") 

Il MSE si appiattisce già per un numero di componenti principali pari a 5.

coefplot(pcr_model) #il modello con 7 componenti principali ha un coefficiente di regressione negativo molto alto 

il modello con 7 componenti principali ha un coefficiente di regressione positivo alto

pcr_model <- pcr(as.numeric(Cfidelity )~., ncomp = 7,
                 data = dati2,
                 subset= train,
                 scale = TRUE, 
                 validation = "CV")

summary(pcr_model) 

## Data:    X dimension: 15998 24 
##  Y dimension: 15998 1
## Fit method: svdpc
## Number of components considered: 7
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV          0.4746   0.3032   0.3032   0.3031   0.3031   0.2958   0.2953
## adjCV       0.4746   0.3032   0.3032   0.3031   0.3031   0.2958   0.2953
##        7 comps
## CV      0.2920
## adjCV   0.2919
## 
## TRAINING: % variance explained
##                        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## X                        22.09    28.86    35.09    40.99    46.54    51.45
## as.numeric(Cfidelity)    59.19    59.20    59.20    59.24    61.19    61.31
##                        7 comps
## X                        56.26
## as.numeric(Cfidelity)    62.18

I valori 59.19 59.20 59.20 59.24 61.19 61.31 62.18 rappresentano la variabilità percentuale spiegata per la variabile y

pcr.pred= predict(pcr_model ,x[test,] ,ncomp=7)
mean((pcr.pred-y.test)^2) # MSE 0.08625655

## [1] 0.08625655

table(as.numeric(dati_test$Cfidelity), round(pcr.pred,0)) #

##    
##        1    2
##   1 1981  705
##   2   72 5242

777 casi su 8000 sono malclassificati, quindi il modello in un’ottica di previsione funziona bene con un errore di malclassificazione di 794/8000= 0.099 –>9,9%

#PLS

plsr_model <- plsr(as.numeric(Cfidelity)~., data = dati2 , 
                   subset = train,
                   scale = TRUE, 
                   validation = "CV")




summary(plsr_model)

## Data:    X dimension: 15998 24 
##  Y dimension: 15998 1
## Fit method: kernelpls
## Number of components considered: 24
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV          0.4746   0.2967   0.2846   0.2796   0.2789   0.2785   0.2785
## adjCV       0.4746   0.2967   0.2845   0.2796   0.2789   0.2785   0.2784
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## CV      0.2784   0.2783   0.2783    0.2783    0.2783    0.2782    0.2782
## adjCV   0.2783   0.2783   0.2782    0.2782    0.2782    0.2782    0.2782
##        14 comps  15 comps  16 comps  17 comps  18 comps  19 comps  20 comps
## CV       0.2782    0.2782    0.2782    0.2782    0.2782    0.2782    0.2782
## adjCV    0.2782    0.2782    0.2782    0.2782    0.2782    0.2782    0.2782
##        21 comps  22 comps  23 comps  24 comps
## CV       0.2782    0.2782    0.2782    0.2782
## adjCV    0.2782    0.2782    0.2782    0.2782
## 
## TRAINING: % variance explained
##                        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## X                        22.04    27.03    29.35    33.60    36.67    41.46
## as.numeric(Cfidelity)    60.93    64.15    65.40    65.56    65.65    65.68
##                        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
## X                        44.03    47.06    50.97     54.49     56.55     60.54
## as.numeric(Cfidelity)    65.71    65.73    65.73     65.74     65.74     65.74
##                        13 comps  14 comps  15 comps  16 comps  17 comps
## X                         64.95     68.55     70.78     73.15     75.63
## as.numeric(Cfidelity)     65.74     65.74     65.74     65.74     65.74
##                        18 comps  19 comps  20 comps  21 comps  22 comps
## X                         79.14     80.57     83.70     87.47     91.60
## as.numeric(Cfidelity)     65.74     65.74     65.74     65.74     65.74
##                        23 comps  24 comps
## X                         95.79    100.00
## as.numeric(Cfidelity)     65.74     65.74

I valori 60.93 64.15 65.40 65.56 65.65 rappresentano la percentuale di variabilità spiegata per la variabile y. Rispetto al modello PCR risulta leggermente più alta.

lm2<-lm(as.numeric(dati_train$Cfidelity)~plsr_model$scores[,1:10])
summary(lm2) 

## 
## Call:
## lm(formula = as.numeric(dati_train$Cfidelity) ~ plsr_model$scores[, 
##     1:10])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10447 -0.09221  0.01968  0.16968  1.18481 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.6575197  0.0021969 754.491  < 2e-16 ***
## plsr_model$scores[, 1:10]Comp 1  0.1619042  0.0009602 168.607  < 2e-16 ***
## plsr_model$scores[, 1:10]Comp 2  0.0841267  0.0021710  38.750  < 2e-16 ***
## plsr_model$scores[, 1:10]Comp 3  0.0790540  0.0032728  24.155  < 2e-16 ***
## plsr_model$scores[, 1:10]Comp 4  0.0219615  0.0024833   8.844  < 2e-16 ***
## plsr_model$scores[, 1:10]Comp 5  0.0208429  0.0032570   6.399  1.6e-10 ***
## plsr_model$scores[, 1:10]Comp 6  0.0087906  0.0023858   3.685 0.000230 ***
## plsr_model$scores[, 1:10]Comp 7  0.0154078  0.0040857   3.771 0.000163 ***
## plsr_model$scores[, 1:10]Comp 8  0.0084104  0.0029281   2.872 0.004080 ** 
## plsr_model$scores[, 1:10]Comp 9  0.0041540  0.0027964   1.485 0.137441    
## plsr_model$scores[, 1:10]Comp 10 0.0034422  0.0029972   1.148 0.250791    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2779 on 15987 degrees of freedom
## Multiple R-squared:  0.6574, Adjusted R-squared:  0.6572 
## F-statistic:  3067 on 10 and 15987 DF,  p-value: < 2.2e-16

il valore dell’ Adjusted R-squared: 0.657.

validationplot(plsr_model)

# considerando solo le prime 7 componenti...

plsr_model <- plsr(as.numeric(Cfidelity)~., ncomp = 7,
                   data = dati2,
                   subset= train,
                   scale = TRUE, 
                   validation = "CV")


summary(plsr_model)

## Data:    X dimension: 15998 24 
##  Y dimension: 15998 1
## Fit method: kernelpls
## Number of components considered: 7
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV          0.4746   0.2967   0.2846   0.2797   0.2789   0.2785   0.2785
## adjCV       0.4746   0.2967   0.2846   0.2797   0.2789   0.2785   0.2784
##        7 comps
## CV      0.2784
## adjCV   0.2784
## 
## TRAINING: % variance explained
##                        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## X                        22.04    27.03    29.35    33.60    36.67    41.46
## as.numeric(Cfidelity)    60.93    64.15    65.40    65.56    65.65    65.68
##                        7 comps
## X                        44.03
## as.numeric(Cfidelity)    65.71

la percentuale di varianza in fedelt con 7 componenti 65.71. Risulta più alta della percentuale spiegata col modello finale dal modello PCR a parità di componenti principali (7)

lm2<-lm(as.numeric(dati_train$Cfidelity)~plsr_model$scores[,1:7])
summary(lm2)

## 
## Call:
## lm(formula = as.numeric(dati_train$Cfidelity) ~ plsr_model$scores[, 
##     1:7])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11150 -0.09145  0.02034  0.16783  1.19204 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    1.6575197  0.0021975 754.284  < 2e-16 ***
## plsr_model$scores[, 1:7]Comp 1 0.1619042  0.0009605 168.561  < 2e-16 ***
## plsr_model$scores[, 1:7]Comp 2 0.0841267  0.0021716  38.739  < 2e-16 ***
## plsr_model$scores[, 1:7]Comp 3 0.0790540  0.0032737  24.148  < 2e-16 ***
## plsr_model$scores[, 1:7]Comp 4 0.0219615  0.0024840   8.841  < 2e-16 ***
## plsr_model$scores[, 1:7]Comp 5 0.0208429  0.0032579   6.398 1.62e-10 ***
## plsr_model$scores[, 1:7]Comp 6 0.0087906  0.0023865   3.684 0.000231 ***
## plsr_model$scores[, 1:7]Comp 7 0.0154078  0.0040868   3.770 0.000164 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2779 on 15990 degrees of freedom
## Multiple R-squared:  0.6571, Adjusted R-squared:  0.657 
## F-statistic:  4378 on 7 and 15990 DF,  p-value: < 2.2e-16

Il valore di Adjusted R-squared: 0.657 rimane quasi invariato ma questa volta sono state utilizzate solamente le prime 7 componenti principali

pls.pred=predict(plsr_model,x[test,],ncomp=7)


mean((pls.pred-y.test)^2) # MSE = 0.07939584

## [1] 0.07930442

table(as.numeric(dati_test$Cfidelity), round(pls.pred,0)) # 775 malclassificati su 8000. 775/8000= 0.0968 --> 9,68%

##    
##        1    2
##   1 1986  700
##   2   75 5239

Il MSE = 0.07939584. Ci sono stati 775 malclassificati su 8000. 775/8000= 0.0968 –> 9,68%

Conclusioni

Successivamente all’analisi esplorativa del dataset, si è proceduto alla creazione della variabile Cfidelity, ottenuta dalla variabile originale “Fidelity”. In seguito, si è suddiviso il dataset casualmente per 2/3 in train e per 1/3 test set. Sono stati applicati 3 modelli sul training set con i quali si è prevista la varibile Cfidelity. Il modello Lasso ha classificato erroneamente 778 osservazioni sulle 8000 del test set. Il modello creato con PCR ha classificato erroneamente 777 osservazioni. Infine il modello PLS ha classificato erroneamente 775 osservazioni del test set.